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Low-thrust trajectory design is tightly coupled with spacecraft systems design. In partic- 
ular, the propulsion and power characteristics of a low-thrust spacecraft are major drivers in 
the design of the optimal trajectory. Accurate modeling of the power and propulsion behav- 
ior is essential for meaningful low-thrust trajectory optimization. In this work, we discuss 
new techniques to improve the accuracy of propulsion modeling in low-thrust trajectory op- 
timization while maintaining the smooth derivatives that are necessary for a gradient-based 
optimizer. The resulting model is significantly more realistic than the industry standard and 
performs well inside an optimizer. A variety of deep-space trajectory examples are presented. 


INTRODUCTION 


Solar electric propulsion (SEP) is a means of propelling spacecraft that provides very high specific impulse 
(Isp) at the expense of low thrust and high power requirements. The high J,,, of a SEP system enables 
many classes of mission that are difficult or impossible with chemical propulsion, including rendezvous with 
multiple asteroids, large-scale sample return, and delivery of very large payloads to the outer solar system. 
SEP trajectory design is more complex than chemical propulsion trajectory design and both the modeling and 
optimization techniques used in such designs are the subject of ongoing research. 


SEP trajectory design is distinct from chemical propulsion trajectory design because the thrust levels are 
low and the maneuver durations are long. The impulse maneuver approximation does not hold and instead it 
is necessary to integrate the equations of motion across a maneuver arc as shown below. 
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where 7’ is maximum possible thrust force, ™ is available mass flow rate, P is power available to the thruster 
as a function of r, distance from the sun, and the wu; are components of a throttle control vector. 
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Electric propulsion systems are composed of a thruster, a power processing unit (PPU), and a propellant 
feed system. The propulsion system may be controlled by varying the input power, the beam voltage setting 
on the PPU, and/or a valve in the feed system that controls 7. PPUs and feed systems are typically designed 
with discrete choices of voltage and m in mind. The throttle table for an electric propulsion system therefore 
resembles a grid in voltage and m. For example the NASA’s evolutionary xenon thruster (NEXT) propulsion 
system ||1] is designed to operate at a range of discrete voltages from 275 to 1800 V and a range of discrete 
mass flow rates from 1.85 mg/s to 5.76 mg/s. As currently designed, the NEXT system must be operated at 
points which appear on the voltage-7m grid, each of which has been tested in laboratory conditions. Each set 
point provides a unique thrust, 7’, and specific impulse, J,,,,. Other SEP systems are similar. 


Thruster performance(T’ and m) on power (P) is a major driver on the shape of the optimal trajectory for 
any SEP spacecraft whose distance from the sun, and therefore available power, change over the course of a 
mission. The more realistic the models for T’ and mm, the more closely the trajectory will resemble one that is 
actually physically acheivable. However there is a conflict between the need for realism in the model and the 
need for a model that is compatible with trajectory design software. 


It is very challenging to model a propulsion system with discrete settings in the context of the gradient- 
based optimizers which are typically used to design low-thrust missions. Both direct and indirect optimization 
methods require a continuously differentiable force model, which conflicts with the discrete nature of the 
actual propulsion systems. Traditionally, low-thrust trajectory optimization software packages reduce the 
2-dimensional throttle grid (each output is a function of both input power and voltage) to a 1-dimensional 
path (each output is a function of input power only). The models typically curve-fit the chosen subset, with 
a well-behaved “smooth” function such as an interpolating polynomial that is conducive to quasi-Newton, 
gradient-based NLP methods. Interpolating polynomials attempt to capture the aggregate behavior of the 
desired subset over the duration of the EP thrusting and moreover do not contain sharp changes in the first 
derivative. Such behavior prevents the thruster model from being the culprit of numerical over-sensitivity 
in the NLP search direction. In first-order methods, such numerical sensitivities can cause convergence 
“chatter” or outright limited robustness due to a lack of resolution in the quasi-Newton step. Smooth models 
are therefore the preferred choice for preliminary design searches, when the thruster firing history is unknown. 
A typical choice is a fifth-degree polynomial as a function of power available to the PPU: 


mm = emo + CmiPa + Cm2P? + ¢m3P3 + CmaPt + emsP2 (5) 
T = cy + Py + CoP? + 43P? + yg Pt + yp PP (6) 


where the c; and c,, coefficients are used to approximate the thruster’s stepped behavior over the range of 
power levels. 


While inexpensive to evaluate and continuously differentiable, the industry-standard polynomial propul- 
sion models are not very accurate. First, the true throttle grid is two-dimensional so one has to filter out 
throttle levels to select a priori which set of points to include in a 1-dimensional polynomial. In a particular 
mission, it may be beneficial to operate in a high-voltage, low-m regime for part of the trajectory and in a 
low-voltage, high-7m regime for another part. Simplifying to a 1-dimensional set could therefore exclude the 
optimal trajectory. Second, the polynomials do not always accurately match the 1-dimensional sets of throttle 
points as shown in Figures[I]and [2] The particular polynomial models and throttle settings shown here are 
those that were provided for the NEXT propulsion system for the Discovery 14 contest (1). The throttle set- 
tings used in these figures and throughout this work are provided by the NASA Glenn Research Center (2). 


The actual throttle grid is discontinuous and therefore non-differentiable and incompatible with a gradient- 
based optimizer. The industry standard polynomial approximation is continuously differentiable but requires 
an a priori assumption about which path to follow through the throttle grid and, even once a path is chosen, 
is inaccurate relative to the actual capabilities of the propulsion system. It is therefore desirable to develop 
new models which: 


1. more accurately model the EP system performance 
2. are continuously differentiable 
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Figure 1: Polynomial thrust model compared to actual throttle grid 
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Figure 2: Polynomial thrust model compared to actual throttle grid 


3. allow the optimizer to choose any point in the 2-dimensional throttle grid 


In this work, novel propulsion models will be presented which satisfy some or all three of the above criteria. 
These new propulsion models are compatible with gradient-based optimizers but closely track the actual 
throttle grid. Multiple versions will be presented: some that require the analyst to choose a 1-dimensional path 
a priori and one that allows a gradient-based optimizer to choose points anywhere in the full 2-dimensional 
throttle grid. The new methods will be demonstrated on realistic mission designs and the results will be 
compared to those found by solving the same problem with the industry-standard polynomial propulsion 
model. 


NON-DOMINATED SET 


First, in order to choose 1-D paths through any arbitrary throttle box, we briefly introduce the concept of 
a non-dominated set. Similar to the formation of a pareto-front in multi-objective optimization, the distinct 
throttle settings of an electric thruster can be filtered to remove any throttle point which is worse than at 
least one other throttle point in both of the filtering objectives. Because all 1-D paths will become functions 
of input power, one of the filtering objectives is always input power. The second filtering objective is used 
to describe the 1-D path chosen. For example, to form the high-m non-dominated set, all throttle settings 
are removed which require both more input power and a higher mass flow rate of propellant than any other 
individual throttle setting. This can also be done to form the low-m set, the low-J,,, set and the high-thrust 
set (the high thrust set is typically similar or identical to the high-m set). These four non-dominated sets for 
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the NEXT propulsion system are shown in Figures[3] [4] and|5] 


Note that these four sets represent the most logical choices for 1-D paths through a throttle grid, however, 
the 1-D models below could be used for any arbitrary set of throttle points that have unique levels of input 
power. The most useful 1-D paths will have monotonically increasing m and thrust (as input power increases) 
or else they will confuse the gradient-based optimizer and reduce its effectiveness. 


POLYNOMIAL MODELS 
Throttle Point Curve Fits 


There are many ways to generate a polynomial fit to a throttle table. The most common approach is 
to perform a least-squares curve fit to a selected set of discrete throttle points. In this case, the predicted 
performance between throttle points is likely over or underestimated. Further, to accurately model the general 
shape of the throttle table, individual throttle points will likely be “missed” by the smooth curve. 


The NEXT engine curve fit provided with the 2014 NASA Discovery Proposal is an example of a smooth 
model that approximated discrete set points (as shown in Figures |1} and [2). However, for trajectories that 
require many thruster set points, this choice results in an over-approximation of the thruster performance 
between steps (see Figure (6). A performance degradation when switching between smooth and more realistic 
stepped models can therefore be expected. 
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Figure 6: NEXT throttle table 11 and 1-D “high-thrust” approximations. 


Further, the Discovery 14 polynomial models miss the throttle points in many locations. This is particularly 
true at low power levels, below 2 kW. This is significant because as a spacecraft moves farther out in the 
solar system, available power drops with 1/r? and so, assuming a reasonable array size, any spacecraft 
going to the main asteroid belt or beyond will spend most of its time thrusting at low power levels. The 
low-power inaccuracy of the “standard” polynomial model will have negative consequences for outer-solar 
system mission design. Of course, one could construct polynomials that are more accurate in the low-power 
regime, but this would likely be at the expense of accuracy in other regions of the throttle grid, and in fact, any 
polynomial fit is likely to be inaccurate for certain ranges of input power. More importantly, any polynomial 
model also allows the propulsion system to operate at settings that do not appear in the throttle grid and are 
not physically realizable. This can result in trajectories that are overly optimistic about performance or worse, 
are also not physically realizable. 


Integral Performance Matching 


Here, we seek to derive a better approach for fitting the available thruster points. We evaluate the high- 
thrust and high-/,,, sets of NEXT throttle levels. Rather than derive a polynomial that matches the throttle 
points identically, we propose a method to approximate the overall behavior. In this approach, any single 
discrete point on the polynomial will have a nontrivial error, but the behavior overall (over many power 
levels) will more accurately model the throttle table. 


A useful polynomial would have the following properties: 


1. The integral or “area under the curve” over the range of available power is roughly preserved. 
2. As possible, the integral over each throttle setting is preserved. 

3. The initial and final throttle points are preserved. 

4. The polynomial is monotonically increasing. 


By preserving the integral of the throttle table, we ensure that the net performance neither over or un- 
derestimates the true performance. That is, the polynomial’s overestimates must cancel its underestimates. 
Ideally, these errors will also be balanced over each throttle setting. Although any single throttle point may 
not be perfectly modeled, it’s useful to capture the behavior at the minimum and maximum power settings. 
The minimum power setting is often encountered in deep space missions, and the maximum power setting is 
often the most efficient and encountered when power is sufficiently high. Any error in these two values could 
have a significant impact on overall accuracy. Finally, it’s favorable that the polynomial be monotonically in- 
creasing because this better emulates the actual thruster behavior, in which increasing power always increases 
the available thrust and specific impulse. If the polynomial has a negative slope, an increase in power would 
result in poorer engine performance. This behavior could cause a trajectory optimizer to incorrectly favor a 
lower available power. 


To this end, we construct a weighted least squares problem to incorporate these behaviors. A polynomial 
fit point is assigned to the middle of each vertical and horizontal step in the throttle table. These can be seen 
in Figure [7] These points are assigned a weight such that their sum is 1.0. An equation for the integral of 
the polynomial is used to constrain the sum of the errors. Two equations are added to constrain the initial 
and final points. Finally, the result is tuned by selectively constraining the slope of the polynomial at key 
points. This manual step prevents large overshoot behaviors. Table [I] gives the sets of least squares weights 
and relevant parameters. This table also includes a description of the constraint equations used. 


The results of this process are given in Figures [710]and Tables [2] and [3] The performance of these poly- 
nomials is related to the polynomials presented in the NASA Discovery 2014 AO document. The custom 
polynomials derived here better balance the errors, where the NASA Discovery AO fits typically overestimate 
performance. The utility of these polynomials is further evaluated in two trajectory optimization examples 
given in sections below. 


Table 1: Polynomial Weights and Inputs (based on EOL data) 


High Thrust High Isp 

Parameter Thrust Mass Flow Rate Thrust Mass Flow Rate 

Throttle Levels Incorporated —_1-10,13,14,18,19,23,24,28,29,33,37-40 1-12,17,22,40 

Polynomial Order Selected 4 6 6 6 

Number of Interior Fit Points 44 16 28 8 

Interior Fit Point Weights 1.0/44 1.0/16 1.0/28 1.0/8 

End Point Weights 25.0 25.0 25.0 25.0 

Integral Equality Weight 0.01 0.01 0.01 0.01 

Additional Slope Weight None 10.0 10.0 10.0 

Additional Slope Constraints 3.0E-07 kg/s/kW @ 0.653 kW | 0.005 N/kW @ 5.829 kW _ 1.0E-09 kg/s/kW @ 0.653 kW 

1.0E-09 kg/s/kW @ 7.360 kW 1.0E-08 kg/s/kW @ 0.806 kW 

2.0E-07 kg/s/kW @ 2.132 kW 
2.0E-07 kg/s/kW @ 5.821 kW 
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Figure 7: Polynomial Fits to the NEXT Thrust Throttle Levels 
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Figure 8: Percent Error in Polynomial Fits to NEXT Thrust Levels 
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Figure 9: Polynomial Fits to the NEXT Mass Flow Rate Throttle Levels 
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Figure 10: Percent Error in Polynomial Fits to NEXT Mass-Flow Rate Levels 


Table 2: Polynomial Fits to the NEXT Thrust Throttle Steps (based on EOL data) 


High Thrust High Isp 
Parameter NASA D14 AO Custom NASA D14 AO Custom 
cro (N) 1.19388817E-02 1.53147067E-02 3.68945763E-03 9.06756030E-03 
C41 (N/kW) 1.60989424E-02 7.77083385E-03 4.054325 10E-02 


C2 (N/kW?) 
13 (N/ kw?) 
cra (N/kW*) 
cis (N/kW?®) 


1.14181412E-02 
-2.04053417E-03 
1.19388817E-02 


1.30817266E-02 
-2.22362321E-03 
1.15864108E-04 


-7.91621814E-03 
1.72548416E-03 
-1.11563126E-04 


7.96175232E-03 
4.22024750E-02 
-3.02627265E-02 
9.285 14190E-03 


0.0 0.0 0.0 -1.29785608E-03 
cg (N/KW®) 0.0 0.0 0.0 6.75040037E-05 
Integral Error +6.24% +0.34% +28.87% 4.25% 
Initial Point Error +5.51% +0.00% +7.39% 0.00% 

Final Point Error -0.13% -0.00% -0.32% 0.00% 
Max Overestimate = 28.2% (@2.14kW) = 15.1% (@2.14kW) | 78.4% (@7.35kKW) 77.7% (@7.35kW) 
Max Underestimate 4.5% (@0.94kW) = 12.9% (@0.94kW) | 5.4% (@1.30kW) = 15.3% (@4.31kW) 


Table 3: Polynomial Fits to the NEXT Mass-Flow Rate Throttle Steps (based on EOL data) 


High-Thrust High-Isp 
Parameter NASA D14 AO Custom NASA D14 AO Custom 
Cmo (kg/s) 2.75956482E-06 1.39733478E-06 2.22052 155E-06 


Cmi1 (kg/s/kW) 

Cm2 (kg/s/kW?) 
Cm3 (kg/ s/kW*) 
Cma (ke/s/kW*) 


-1.71102132E-06 
1.21670237E-06 

-2.07253445E-07 
1.10213671E-08 


1.83615417E-06 
-2.08514759E-06 
1.13690887E-06 
-2.61967933E-07 


-1.80919262E-07 
2.77715756E-08 
2.98873982E-08 
-2.91399146E-09 


2.00382111E-06 
-1.06921486E-07 
1.74340894E-07 
- 1.433045 17E-07 
6.65953023E-08 


Cms (kg/s/kW°) 0.0 2.72340756E-08 0.0 -1.22553728E-08 
cme (kg/s/kW®) 0.0 -1.06262899E-09 0.0 7.64533395E-10 
Integral Error +4.81% -0.68% +24.42% 4.43% 
Initial Point Error +6.39% 0.00% +7.22% 0.00% 

Final Point Error -0.05% 0.00% -0.52% 0.00% 

Max Overestimate 43.3% (@2.20kW) 20.2% (@2.21kW) | 73.2% (@7.35kW) 73.0% (@7.35kW) 


Max Underestimate 


2.4% (@5.04kW) 


11.0% (@3.07kW) 


3.8% (@1.22kW) 


10.8% (@4.31kW) 
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Figure 11: Heaviside step function and approximations 


STEPPED MODELING 
1-D Stepped Model 


This section outlines a method to create a mathematical equation which has discontinuous increases from 
one throttle setting to another. This is possible using the so-called Heaviside step function (see equation [7). 
A Heaviside function has a value of one for an input greater than or equal to a critical value, and zero for 
all other inputs. Each throttle point can be modelled with a Heaviside function that turns its value of thrust 
and mass flow rate “on” at its input thruster power. However, the magnitude of each Heaviside function only 
reflects a step increase over the previous throttle point’s value. Let P, be the continuously defined power 


available to the PPU, and x; € {21,72,...,@n} denote either a mass flow rate or thrust magnitude set point 
associated with a power level P; € {P,, P2,..., Py}. 
1 if S$; >0 
=f 4 if S; <0 7) 
5S;=P,—P,; (8) 


The mass flow rate or thrust magnitude available at a given power available is then: 
c= 2A + Soa; = Mena) Hi (9) 
i=2 


Implementing equation[9]solves the first of the three criteria set out above, in that the model output perfectly 
matches the capability of the thruster at any input power. However, this model does not meet the second 
criteria. The differentiability problem is solved by replacing the discontinuous Heaviside functions above 
with a logistics function (Figure[I1) which smoothly and continuously approximates step increases: 


— 1 


i — > oe ! 
1+ 67k: 


The parameter k, or throttle sharpness, determines how closely the logistics function models the step increase 
of a Heaviside function. 


The net model which solves the first two criteria set out above we name the 1-D stepped model and is 
presented in equation : 


w= 0A, + So (a; - 2-1) Hi (11) 
i=2 


10 


and the relevant analytical derivative is: 


dx 5 _hg. ” —~ > _ks, 
ap, = Reilite orl Apnea ye Nes (12) 


For a conceptual example of how this model works, let throttle levels 1 and 2 are consecutive members 
of a chosen non-dominated set requiring | kW and 2 kW to operate and generating 100 mN and 150 mN 
of thrust, respectively. As the power available to the thrusters increases past 2 kW, the Heaviside function 
for throttle level 2 switches “on”. However, because the engine is already outputing 100 mN from throttle 
point 1, the step increase associated with throttle level 2’s Heaviside function is only 50 mN resulting in a net 
output of 150 mN. The input power has not yet been reached for any of the other throttle levels, so they do 
not contribute to the “net” thrust and mass flow rate, even though they are summed over. 


For values of k greater than 1000, the Heaviside step function and the logistics function approximation 
are nearly indistinguishable. However, this has negative effects for the differentiability of the function. As k 
increases, the derivatives with respect to power in the “flat” regions grow closer and closer to zero, and the 
derivative for P, = P; becomes larger and larger. In practice, the value of & should be chosen for the specific 
throttle box, and the capability of the gradient-based optimizer. Further, there could be a separate value of k 
chosen for each throttle level, if so desired, but here we assume that one identical value is used everywhere. 


1-D Stepped Model with Smoothing Radius 


One drawback to the 1-D method described above is that the parameter to control the steepness of the 
discontinuity approximation, k, is non-dimensional. It is often preferable to have this parameter in units of 
the independent variable, input power. This can be accomplished with a different Heaviside approximation. 
If we define a smoothing radius, A, then the Heaviside approximation becomes: 


7 Fin 
A; = ——— 13 
Fyit+ Fi ce 
Fy 4 acs o( sare) S; > —A/2 (14) 
0 S; < —A/2 
fS= (ars) S; < A/2 (15) 
0 S; > A/2 


For gradient-based algorithms, the derivative of equation with H replacing H is readily available as: 


n 


dx dH, dH; 
SOE ag A) =e 1 
dP, dP, pa ti-1) op, ue) 
OF i OF ;,2 
OH; 7 ( as )Fi2 — ( a )Fi1 (17) 
OP (Fi + Fi2)? 
: 1 sas 
OF 4 = cathe FA ) S; > —A/2 (18) 
oP 0 5; < -A/2 
OF ,2 e aitsgrel a) Si < A/2 (19) 
oP 0 S;>A/2 


In contrast to the sharpness, k, approach, this approximation allows direct control over the power smoothing 
radius A of the set points. Choice of an appropriate value is again a compromise between accuracy and 
numerical sensitivity. A value of A = 0.1 kW in trial problems generally satisfies these competing objectives. 
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2-D Stepped Model 


The 1-D model solves the first two of our objectives, but not the third, that of allowing the optimizer to 
choose any throttle level within the full grid. This will first necessitate the introduction of an additional 
control variable. Any parameter which has a discrete value at each throttle setting could be used (thrust, 
Isp, voltage, current, 7), but the fewer settings that the control variable has, the more effective it will be. 
In this work, we will continue to use the NEXT thruster as an example, and as such, 7m, voltage, or current 
are the most logical control variables. As any of these control variables could be used (with slightly varying 


effectiveness), we will simply use Ucommana to represent whichever control variable was selected. 


Similiar to the 1-D model, the 2-D model will use the logistics function approximation to the Heaviside 
step function to turn each throttle point “on”, however in the 2-D model, each throttle point will use Heaviside 
functions to turn the throttle level “off” as well. Because of the increased dimensionality, there is no a priori 
order to the throttle points. From a given pair of P, and Ucommand, depending which (or both) control 
variable changes to move in the throttle grid, many different throttle levels could be the next to be activated. 


Again, the net thrust output by the EP system is a summation of the thrust and mass flow rate functions for 
all throttle points. Assuming that the command variable is voltage, the thrust and mass flow rate generated is: 


t= a «i(Hp,,,., _ Fe og ves _ Ties) (20) 
t=1 


Recalling that Hp was defined in equation|10 we simply now add a more detailed subscript to distinguish 
between the command variable and power available. Hy takes the same form, simply with critical values of 
voltage instead of power: 

— 1 


fy, ~ 1+ e7k(Uvoltage —Vi) 


(21) 


As the throttle table is initally parsed, the values of Poni, Pos fis Moni, and Moff; are found for each 
throttle level. The “on” value for each throttle level always corresponds to the operating conditions of the 
given throttle level. Power and the command variable are sorted in slightly different manners to determine 
the “off” condition. To explain, continue assuming that voltage is used as the command variable. All of the 
throttle points are sorted by their required voltage input. For a given throttle level, the critical value of voltage 
which turns it “off” by way of a negative Heaviside step function is the next greater voltage level, regardless 
of required input power. On the other hand, the critical available power is found by sorting only those throttle 
levels with identical input voltage and selecting the next greater input power. Given both sorting methods, 
the 2-D model will be most effective by selecting a command variable with the fewest discrete levels. For 
example, even though the NEXT thruster has 40 total throttle settings, there are only 12 potential values of 
voltage and 8 potential values of mass flow rate. 


Finally, the derivatives with respect to decision variables are defined in the following equations (again with 
voltage as the command variable): 


Ox “ Fy A —k(Pa—Pon,i) F72 —k(Pa—Poff,i) F72 
aP, = bin. > Hy,5,;)(€ \ , 1D — eM oe HP, ,) (22) 
0 ” = — : = Na 
an = kS > «i(Ap,,.; = FL, ,)(e P(tvorrae— ont) FF? - eB Uvottage— oft AY) (23) 
Uvoltage j=l : , 
CASE STUDIES 


Ceres 


As a first example, consider a hypothetical 2024 Earth-Mars-Ceres rendezvous with the assumptions 
of Table The results of the trajectory optimizations are qualitatively very similar, however the over- 
approximation of the standard polynomial model creates a major mismatch in trajectory performance against 
all alternative propulsion models (see Figure[I2]and Table[5). 
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Table 4: Assumptions for Earth-Mars-Ceres Rendezvous Mission 


Option Value 

Optimization Objective max. final mass 

Launch date boundaries 1/31/2024 — 10/30/2024 
Flyby date boundaries 7/30/2025 — 9/28/2025 
Flight time upper bound 7 years (2555 days) 
Propellant upper bound 500 kg 

Arrival condition at Mars flyby 

Arrival condition at Ceres rendezvous 

Launch vehicle Atlas V 401 

Launch asymptote declination bounds [—28.5°, 28.5°] (Kennedy Space Center) 
Post-launch coast duration 30 days 

Pre-flyby coast duration 30 days 

Solar array Po 15 kW 

Solar array coefficients +; [1.0, 0, 0, 0, 0] 
Spacecraft power coefficients as/.—Cs/- _[1.0, 0, 0] 

Propulsion system varies by trial 

Duty cycle 90% 

Thrust scale factor 100% 

Power margin 0% 

Number of time steps per phase 100 


All trials nearly reach the propellant throughput limit, but the stepped-models (Trials #2 and #3) are less 
capable, leading to a lower launch mass to push the spacecraft through a similar trajectory. The over esti- 
mation is clear from Figure as the stepped models show that it is impossible to fully utilize all of the 
available power. The mismatch becomes most significant during the period following the Mars flyby, where 
significant periods of thrusting at lower throttle settings occur. The gradual shifting between throttle points 
exacerbates the modeling disagreement between the smooth and stepped models. 


The 2-D stepped model (Trial #3) provides only roughly a 1 kg improvement over the 1-D stepped model 
in both delivered mass and propellant usage. The throttle settings are somewhat different however, as the 
optimization algorithm chooses some throttle settings from the high-Isp set which are not available for the 
1-D optimization (see Figure [13). The differences can also be seen with careful observation of Figure[14} as 
the 2-D solution thrusts for slightly longer on both phases of the journey and utilizes increased power, but 
due to the higher Isp of the throttle levels used, less propellant is consumed. In this instance, the 2-D solution 
was found by using the 1-D stepped solution as an initial guess, and the initial guess was actually feasible 
for the 2-D optimization run, as the command variable values were explicitly calculated to match the 1-D 
solution. It is therefore not surprising that the 2-D solution matches so closely to the 1-D solution and in 
fact, only one additional throttle setting was used, as seen in Figure [13] In order to retain the utility of the 
initial guess, the throttle sharpness was kept at k = 10000. As mentioned this causes the derivatives to be very 
small in most regions of the design space and very large as the throttle points switch. The net effect is that 
the gradient based NLP solver has numerical difficulty making significant improvement away from the initial 
guess. Further improvement could likely be found with an improved global search capability, however, this 
does not change the slight shortcoming of the 2-D model. As the throttle sharpness parameter, k, decreases 
or the smoothing radius, A, increases, the 1-D stepped model decreases in accuracy at a much slower rate 
than does the 2-D model. Because the 1-D model only uses step increases from the previous throttle point, it 
is much less sensitive to k. This situation could likely be avoided by developing initial guesses with a lower 
value of k, however this would require increased computation time. 


The fourth trial, utilizing the equal area integral approximation of the high-thrust throttle points in the 
stepped model, significantly improves the trajectory agreement with the 1-D stepped model, including the 
critical event dates, launch C3, and most importantly, the optimal final mass which disagrees with the 1-D 
stepped model by only 3.95 kg. 
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Ceres Rendezvous 


Event #1 
Earth Launch 


Event #2) 0 5 ae ee 


Figure 12: Earth-Mars-Ceres trajectory example. Due to qualitative similarities, only trial #1 is shown. 


Trial 1 Trial 2 Trial 3 Trial 4 
2014 Discovery AO 1-D stepped : 
Engine Model Hahnne high-thrust, A = 0.1 sl, ie eee 
é 10000 polynomial 
polynomial kW 
Launch Date 9/26/2024 9/24/2024 9/25/2024 9/25/2024 
Launch C3 10.55 km?/s? 11.09 km?/s? 11.26 km?/s? 10.74 km?/s? 
Launch DLA 6.1° —1.4° —.85° 1.6° 
Launch Mass 2057.84 kg 1976.04 kg 1976.70 kg 1979.99 kg 
Mars Flyby Date 7/30/2025 7/30/2025 7/30/2025 7/30/2025 
Mars Flyby vo. 2.46 km/s 2.70 km/s 2.70 km/s 2.66 km/s 
Mass at Flyby 1993.05 kg 1895.34 kg 1896.47 kg 1895.66 kg 
Mars Flyby Altitude 300 km 300 km 300 km 300 km 
Ceres Arrival Date 9/25/2031 9/23/2031 9/24/2031 9/24/2031 
Total Flight Time 2555 days 2555 days 2555 days 2555 days 
Mass After Rendezvous 1557.84 kg 1476.04 kg 1477.38 kg 1479.99 kg 
Total Propellant Use 500.00 kg 500.00 kg 499.32 kg 500.00 kg 
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Table 5: Trajectory data for the Earth-Mars-Ceres example 


Thrust, % 


ONEXT X2D Stepped Model 


41D Stepped Model ™ 2014 Discovery AO High Thrust Polynomial 
© Equal Area Integral Polynomial 
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Figure 13: Comparison of engine model outputs for the sample Ceres mission. Nondimensional thrust is 
plotted, scaled by the value of the AO polynomial at any given mass flow rate. 
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Figure 14: Comparison of power usage for the four trials 
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Trojan Tour 


A notional SEP mission to the Trojan asteroids Nauplius and Odysseus is used as the second benchmark 
case for the new propulsion models. This benchmark is particularly timely because the “Trojan Tour and 
Rendezvous” is listed as one of the acceptable targets in the New Frontiers 4 Announcement of Opportunity 
that was recently released by NASA (3). As a mission of current interest, the Trojan Tour and Rendezvous 
is an ideal case with which to validate a new design tool that may be used to propose such missions. The 
Announcement of Opportunity states that the Trojan Tour and Rendezvous “is intended to examine two 
or more small bodies sharing the orbit of Jupiter, including one or more flybys followed by an extended 
rendezvous with a Trojan object.” Accordingly, our benchmark mission performs a flyby of the magnitude 
10.7 Trojan 9712 Nauplius before rendezvousing with 1143 Odysseus, which at magnitude 7.93 is one of the 
largest Trojans. This benchmark is appropriate for comparing the NEXT “stepped” propulsion models to the 
industry-standard polynomials because it exercises the propulsion model at a wide range of available power. 


The assumptions for the benchmark mission to Nauplius and Odysseus are listed in Table [6] The sec- 
ond benchmark mission was designed using the Evolutionary Mission Trajectory Generator (EMTG) (41/5). 
EMTG uses a form of global search known as monotonic basin hopping to generate initial guesses, but each 
individual optimization run requires a gradient based NLP solver. For this reason, it is crucial that changes in 
the engine model performance be smooth with continuous derivatives. The mission was designed using both 
polynomial models and both the 1-D and 2-D stepped models. 


Table [7] lists the major events in all four versions of the mission, and Figure [15] is a plot of the optimal 
trajectory. Note that all four optimal trajectories look qualitatively similar in a Universe view, so only Trial 
#2 (1-D stepped model) is shown. Finally, Figure [16] displays the EP system outputs used in each of the 
trajectories. 


As expected, the solution found using the industry standard polynomial severely overestimates the dry mass 
that can be delivered to an Odysseus rendezvous by roughly 9%. On the other hand, in this case, the improved 
polynomial model, while much more accurate, still over estimates the 1-D NEXT system capability by 1.8%, 
but the 2-D capability by only .08%. The equal area integral polynomial is clearly a preferable polynomial 
model for early stage design. 


The 2-D stepped model allows more mass to be delivered to the Trojan rendezvous than the 1-D stepped 
model, and the solutions are 1.7% separated. An improvement is to be expected, as the optimizer has an 
increased number of throttle settings to choose from. However, it is difficult to predict the amount of im- 
provement. Unlike the Ceres mission, the optimizer selected many throttle settings that were not on the 1-D 
non-dominated path. This is primarily due to the spacecraft having two thrusters onboard. Given that the 
solar arrays were sized in order to allow maneuvering at 5 AU, near the Trojan asteroid cloud, there is ex- 
cess power available to fire two thrusters when the spacecraft is still near Earth, and with two thrusters, the 
spacecraft is able to get sufficient thrust capability at the high Isp settings, rather than needing to use solely 
the high thrust settings. 


One unexpected result from this comparison study is that trials 2,3 and 4 all found optimal solutions which 
underloaded the Atlas 551 launch vehicle by roughly 2000 kg. Further, despite launching with lower C3, both 
stepped solutions require less propellant than the trial #1 solution. This is not necessarily surprising for the 
2-D solution which has access to a greater number of throttle control settings, some of which are more fuel 
efficient, but this is unexpected for the 1-D solution. In many cases, switching from a polynomial solution to 
a stepped solution not only decreases the mass deliverable to the final orbit, but also increases the propellant 
required. 


While the solution for all four trials activated the total flight time constraint, the date of Nauplius intercept 
varies by over one month and demonstrates the greatest variability in the four trajectory timelines. Given 
the variability, yet the qualitative similarity in the trajectories, the overall trajectory does not appear to be 
highly sensitive to the flyby date. A further constraint could likely force the timelines to match exactly, with 
a similar variation in mass history. 
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Table 6: Assumptions for the benchmark mission to Nauplius and Odysseus 


Option Value 

Launch window open date 1/1/2024 
Launch window close date 12/31/2024 
Flight time upper bound 12 years 

Arrival condition at Nauplius flyby 

Arrival condition at Odysseus rendezvous 
Launch vehicle Atlas V 551 
Launch asymptote declination bounds [—28.5, 28.5] (Kennedy Space Center) 
Post-launch coast duration 60 days 
Pre-flyby coast duration 30 days 
Post-flyby coast duration 30 days 

Solar array PO 40 kW 

Solar array coefficients 7; [1, 0,0, 0, 0] 
Spacecraft power coefficients as/¢—Cs/c [0.8, 0, 0] 
Propulsion system 2 NEXT engines 


Throttle Logic 

Duty cycle 

Thrust scale factor 

Power margin 

Number of time steps per phase 


minimum number of thrusters 
90% 

99% 

15% 

20 


Event # 3: 
LT rendezvous 
Odysseus 


Event # 2: 
intercept 
Nauplius 


Event # 1: 
launch 
Earth 


Figure 15: Optimal Trojan Tour Mission (all 4 trials qualitatively similar) 
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Trial 1 Trial 2 Trial 3 Trial 4 


Engine Model sae sa 1-D stepped 2-D stepped, k = equal area integral 
: high-thrust, k = 1000 10000 polynomial 
polynomial 
Launch Date 9/15/2024 10/9/2024 9/26/2024 8/10/2024 
Launch C3 12.29 km?/s? 9.79 km?/s? 10.38 km?/s? 8.16 km?/s? 
Launch Mass 3107.91 kg 2935.94 kg 2954.63 kg 2999.79 kg 
Nauplius Intercept Date 5/12/2033 6/15/2033 5/26/2033 5/1/2033 
Nauplius Flyby v.. 2.19 km/s 2.21 km/s 2.20 km/s 2.20 km/s 
Mass at Intercept 1981.59 kg 1831.85 kg 1861.15 kg 1853.23 kg 
Odysseus Arrival Date 9/15/2036 10/9/2036 9/26/2036 8/10/2036 
Total Flight Time 4383 days 4383 days 4383 days 4383 days 
Mass After Rendezvous 1789.21 kg 1641.01 kg 1669.70 kg 1671.17 kg 
Total Propellant Use 1318.70 kg 1294.94 kg 1284.93 kg 1328.62 kg 


Table 7: Trajectory data for the sample Trojan tour mission 


ONEXT X2D Stepped Model 
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Figure 16: Comparison of Engine Model Outputs for the Sample Trojan Tour Mission. Nondimensional 
thrust is plotted, scaled by the value of the AO polynomial at any given mass flow rate. 
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CONCLUSIONS 


Multiple models have been presented for improved modeling of low thrust engine performance for use with 
gradient based optimizers. The standard method for modeling engine performance is to create a polynomial 
for thrust generated and mass-flow rate of the engine as a function of input power only. This can over or 
under estimate performance, predict thruster outputs which are physically unrealistic, and often do not even 
match the throttle points that they are approximating. 


For preliminary design, in order to achieve robust optimization capability, polynomials are useful nonethe- 
less. However, better approximation of the overall performance would be acheived with a modified polyno- 
mial, such as the area integral matching method. Further, using step function approximations, optimization 
is possible using either a 1-D non-dominated sets of the throttle table or the full 2-dimensional table itself. 
The numerical challenges of generating the 2-D set however, makes it much more sensitive to the throttle 
sharpness, and good initial guesses are still required. 


Examples were shown of realistic mission design using each of the models. Polynomial models do not 
accurately model performance, but improved models can get closer than the industry standard least squares 
regression. Further the 2-D stepped model shows improvement over the 1-D throttle sets, as the optimizer is 
given increased control authority. 


Further testing of each stepped model is planned in order to quantify convergence and demonstrate them 
on a wider set of problems. 
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